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Executive  Summary 


Given  measured  water-leaving  radiances  or  remote-sensing  reflectances,  it  is  not  possible  to 
invert  the  radiative  transfer  equation  (RTE)  to  obtain  the  water  column  absorption  and  backscatter 
coefficients  unless  additional  assumptions  are  made,  or  additional  information  is  added,  to  constrain 
the  inversion.  One  approach  to  constraining  the  inversion  is  to  model  selected  unknowns — ^the  so 
called  shape  factors  and  related  quantities — ^as  well  as  possible  in  terms  of  known  quantities,  and 
then  treat  the  shape  factors  as  known  in  the  inversion  process.  This  reduces  the  number  of  unknowns 
in  the  inversion,  at  the  expense  of  introducing  uncertainties  via  the  imperfectly  modeled  shape 
factors.  Clearly,  the  ultimate  success  of  the  inversion  then  depends  on  how  well  the  shape  factors 
can  be  modeled.  This  paper  presents  preliminary  models  for  several  shape  factors  and  related 
quantities  that  are  needed  in  the  inversion  of  the  radiative  transfer  equation.  The  models  presented 
here  were  developed  using  a  data  base  of  Hydrolight-generated  synthetic  data  for  homogeneous  Case 
1  waters  whose  inherent  optical  properties  (lOPs)  are  described  by  simple  bio-optical  models  for  the 
absorption  and  scattering  coefficients. 

If  lOPs  such  as  the  absorption  and  backscatter  coefficients  are  recovered  as  part  of  an  iterative 
inversion  algorithm,  then  the  recovered  lOPs  can  be  used  to  model  the  shape  factors  because  both 
the  recovered  lOPs  and  the  lOP-dependent  shape  factors  are  iteratively  improved  during  the  course 
of  the  RTE  inversion.  In  the  present  modeling,  we  considered  the  absorption  coefficient  a,  the  beam 
attenuation  coefficient  c,  and  the  backscatter  coefficient  6*  to  be  the  recovered  (i.e.,  known)  lOPs. 
The  newly  developed  lOP-dependent  shape  factor  models  can  predict  the  shape  factors  with 
accuracies  ranging  typically  from  two  to  20%.  Using  the  modeled  shape  factors,  we  are  able  to 
predict  the  in-air  remotely  sensed  reflectance  (RSR,  which  is  equivalent  to  the  remote  sensing 
reflectance  RJ  to  within  20%  of  the  correct  (Hydrolight  computed)  values  96%  of  the  time  for  the 
synthetic  data  used  to  determine  the  shape  factor  models,  and  to  predict  RSR  to  within  ±0.0005  sr'* 
86%  of  the  time.  These  results  are  very  encouraging  and  warrant  the  continued  development  of  the 
RTE  inversion  algorithm  and  its  evaluation  using  both  synthetic  (Hydrolight  generated)  data  and  real 
data  obtained  at  sea. 

The  details  of  this  work  are  presented  here  in  a  format  suitable  for  submission  to  the  Journal  of 
Geophysical  Research  as  the  first  part  of  a  two-part  paper.  The  present  Part  I  presents  the  underlying 
theory  and  models  for  the  shape  factors  and  related  quantities.  Part  n  will  present  the  results  of  the 
application  of  the  RTE  inversion  algorithm  to  both  synthetic  and  real  data. 
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ABSTRACT 


We  first  show  that  the  in-water,  shape-factor  formulation  of  the  radiative  transfer  equation  (RTE) 
can  be  modified  to  yield  exact  in-air  expressions  for  the  remote  sensing  reflectance,  /?„,  and  the 
equivalent  remotely  sensed  reflectance,  RSR^.  Our  form  of  the  RTE  can  be  configured  for  inherent 
optical  property  retrievals  via  standard  linear  matrix  inversion  methods.  Linear  matrix  inversion  of 
the  shape-factor  RTE  is  exact  in  the  sense  that  no  approximations  are  made  to  the  RTE.  Thus  errors 
in  retrieved  inherent  optical  properties  (lOPs)  are  produced  only  by  uncertainties  within  (1)  the 
models  for  the  shape  factors  and  related  quantities  and  (2)  the  lOP  models  required  for  inversion. 
We  then  use  Hydrolight  radiative  transfer  calculations  to  derive  analytical  models  for  the 
backscattering  shape  factor,  radiance  shape  factor,  fi-actional  forward  scattering  coefficient,  ratio  of 
air-to-water  mean  cosines,  and  diffuse  attenuation  coefficient  for  in-water  upwelling  radiance  for 
Case  1  waters.  These  shape  factor  models  can  predict  the  various  shape  factors  with  accuracies 
ranging  typically  fi-om  two  to  20%.  Using  the  modeled  shape  factors,  we  are  able  to  predict  the  in-air 
remotely  sensed  reflectance  RSR^  to  within  20%  of  the  correct  (Hydrolight-computed)  values  96% 
of  the  time  for  the  synthetic  data  used  to  determine  the  shape  factor  models,  and  to  predict  RSR^  to 
within  ±0.0005  sr '  86%  of  the  time.  Error  propagation  matrices  suggest  that  inherent  optical 
property  retrievals  are  primarily  influenced  by  uncertainties  in  the  backscattering  shape  factor,  rather 
than  by  imcertainties  in  the  remaining  shape  factors  and  their  related  components.  Our  model  for 
the  backscattering  shape  factor  gives  predictions  that  are  correct  to  within  5%  for  96%  of  the 
synthetic  data  points  used  to  develop  the  model. 


1.  Introduction 


Semianalytic  radiance  models  [Gordon  et  al.,  1988;  Morel  and  Gentili,  1996]  can  be  readily 
inverted  by  linear  matrix  methods  [Hoge  et  al.,  1999a,b]  to  provide  oceanic  inherent  optical 
properties  (lOPs).  Such  inversions  are  well  conditioned  [Hoge  and  Lyon,  1996]  and  promise  a 
powerful  method  of  simultaneously  retrieving  constituent  absorption  and  backscattering  coefficients 
in  the  upper  surface  layer  of  the  world's  oceans  using  satellite  data.  However,  semianalytic  radiance 
models  do  not  provide  an  exact  framework  to  account  for  all  possible  environmental  and  viewing 
conditions  [Weidemann  et  al.,  1995]. 

The  radiative  transfer  equation  (RTE)  can  provide  exact  inverse  solutions,  but  the  RTE  is  not 
easily  inverted  for  many  remote  sensing  situations  [Zaneveld,  1995].  We  therefore  investigate  the 
inversion  of  a  specific  form  of  the  RTE,  namely  a  modified  version  of  the  shape-factor  formulation 
of  Zaneveld  [1995].  Some  of  the  motivation  for  this  work  comes  from  the  distinct  need  for  highly 
accurate  methods  to  retrieve  the  absorption  coefficients  of  the  chlorophyll  accessory  pigment 
phycoerythrin  [Hoge  et  al.,  1999b].  To  this  end  the  absorption  coefficients  of  chlorophyll  and 
CDOM  must  be  accurately  retrieved,  otherwise  weaker  absorbing  constituents  (such  as 
phycoerythrin)  will  be  obscured. 

In  this  paper  we  (1)  show  that  the  shape-factor  form  of  the  RTE  can  be  readily  configured  into 
linear  form  for  simultaneous  retrieval  of  oceanic  lOPs  using  standard  matrix  methods;  (2)  derive  the 
RTE  inversion  for  the  principal  "big  three"  lOPs,  namely  the  phytoplankton  absorption  coefficient, 
the  CDOM+detritus  absorption  coefficient,  and  the  total  constituent  backscattering  coefficient;  (3) 
develop  shape  factor  and  related  models  required  for  inversion,  namely  backscattering  and  radiance 
shape  factors,  the  diffuse  attenuation  coefficient  for  upwelling  radiance,  the  ratio  of  average  cosines 
of  the  air  and  water  downwelling  irradiances,  and  the  fractional  forward  scattering  coefficient;  and 
(4)  assess  the  propagation  of  errors  into  the  lOP  state  vector  resulting  from  errors  in  the  data-model 
matrix  and  hydrospheric  vector  as  well  as  shape  factor  and  related  models.  Our  ultimate  objective 
is  to  determine  if  the  shape-factor  RTE  matrix  inversion  methodology  will  result  in  accurate  and 
efficient  algorithms  for  application  to  satellite  ocean  color  data.  The  present  paper  presents  the 
underlying  theory  and  develops  the  needed  models  for  the  shape  factors  and  related  quantities.  The 
companion  paper  [Hoge  et  al,  2002]  applies  this  theory  to  synthetic  and  real  data. 
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2.  The  Shape  Factor  Form  of  the  Radiative  Transfer  Equation 


In  a  plane  parallel  mediiun  without  internal  sources  or  inelastic  scattering,  the  radiative  transfer 
equation  is 


dl(e,(p,z)  _ 


-c(z)I(0,(p,z)  +  f  jfP(0,(p;0V;2)  1,(0 V^)  sine' (1) 


(See  Notation  section  for  definition  of  symbols.)  Zaneveld  [1995, 1982]  showed  that  the  Eq.  (1)  can 
be  rewritten  in  terms  of  the  in-water  remotely  sensed  reflectance  (RSR)  as 


RSR  = 


-cos0^0,(p,z)  +  c(z)  - 


where  the  dimensionless  backscatter  shape  factor  is  given  by 


j  j  p(0,(p;e',q)',z)  I^0',(p',z)  sin  Q'dB^  dip' 


/^(e,(p,z)  = 


the  dimensionless  radiance  shape  factor  is  given  by 


/,(0,(p,z)  = 


i  j  P(0,(p;0',(p',z)l„(0',(p',z)sin0'i/0'c/(p' 

0  Jt/2 _ 

6/z)Lf(z) 


and  the  diffuse  attenuation  fimction  for  upwelling  radiance  ^0,(p,z),  with  units  of  1/meters,  is  given 
by 

-  1 


=  - 


dz 


Equation  (2)  is  Zaneveld's  [1995]  Eq.  (7)  and  is  exact  because  it  is  simply  a  restatement  of  the  RTE 
(1)  for  upward  directions  using  definitions  (3)-(5).  Subscripts  d  and  u  on  the  radiance  L  explicitly 
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remind  us  that  the  radiance  in  Eq.  (3)  is  downwelling  (0^0^  7t/2),  whereas  the  radiance  in  Eqs.  (4) 
and  (5)  is  upwelling  {idl  <  0  ^  Jt). 

The  numerator  of  the  ^  shape  factor  of  Eq.  (3)  shows  how  much  downwelling  radiance  is 
scattered  upward  into  direction  (0,(p).  The  denominator  is  the  same  quantity  evaluated  for  the  special 
case  of  an  isotropic  volume  scattering  function  (in  which  case  P  =  ItJAn).  Thus  the^^  shape  factor 
is  a  measru:e  of  how  much  the  actual  phase  function  differs  from  a  constant  over  the  backscattering 
directions.  Similarly,  the  numerator  of in  Eq.  (4)  shows  how  much  up  the  upwelling  radiance  is 
forward  scattered  into  direction  (0,(p).  The  denominator  is  the  same  quantity  evaluated  for  the 
special  case  of  an  isotropic  upwelling  radiance  distribution  whose  magnitude  is  L^,  and  for  the 
special  case  of  an  isotropic  volume  scattering  function.  Clearly,  these  shape  factors  depend  both  on 
the  lOPs  (namely  on  the  volume  scattering  function,  in  this  case)  and  on  the  ambient  radiance 
distribution,  as  does  the  diffuse  attenuation  coefficient  of  Eq.  (5).  These  quantities  therefore  are 
unknown  terms  in  Eq.  (2),  if  Eq.  (2)  is  to  be  inverted  to  obtain  the  lOPs  a  and  from  measured 
upwelling  radiances  and  downwelling  irradiances.  The  fact  that  shape  factors  are  unknown  prevents 
the  RTE  (2)  from  being  inverted  unless  further  assumptions  are  made  about  the  values  of  the  shape 
factors.  Modeling  these  unknowns  in  terms  of  known  quantities  is  the  major  focus  of  this  paper. 

Equations  (l)-(5)  are  valid  at  any  depth  within  an  arbitrarily  stratified  water  column.  Our  interest 
is  in  remote  sensing  of  near-surface  water  lOPs.  We  therefore  need  to  relate  the  quantities  in  Eqs. 
(2)-(5),  when  evaluated  just  beneath  the  mean  sea  surface,  to  quantities  in  air  just  above  the  sea 
surface,  which  can  be  deduced  via  in-air  remote  sensing  techniques.  Equation  (2)  can  be  converted 
into  a  form  suitable  for  above-water  remote  sensing  applications  as  follows.  The  n-squared  law  for 
radiance  transmittance  across  a  boundary  between  two  media  [Mobley,  1 994,  Eq.  (4.2 1)]  can  be  used 
to  convert  the  in-water  upwelling  radiance  just  beneath  the  sea  surface,  I,X6,(p,z=0),  to  the  water¬ 
leaving  radiance  just  above  the  sea  surface,  Z,„„(0j,(p): 

2 

T„(0,(p,Z=0)  =  (6) 


Subscript  a  denotes  values  in  air,  just  above  the  mean  sea  surface;  depth  z  =  0  denotes  values  in 
water,  just  beneath  the  sea  surface.  The  in-air  polar  angle  9g  associated  with  I,„„(03,(p)  is  the  refracted 
viewing  angle  above  the  sea  surface  obtained  by  applying  Snell's  Law  to  the  in-water  angle  0.  The 
downwelling  scalar  irradiance  EJ^z)  is  converted  to  the  downwelling  plane  irradiance  E’/z)  via  the 
mean  cosine  of  the  downwelling  radiance,  p.^:  E^(z)=  Ed(z)/p(,.  The  plane  irradiance  just  beneath 
the  sea  surface  can  be  related  to  the  in-air  value  via  [Mobley,  1994,  Eq.  (7.19)] 
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(7) 


£'d(z=0)  =  £'daT  /  (1  +  rR).  Combining  these  results  gives  [Mobley,  1994,  Eq.  (10.27)] 

1 1  Eja 


In  Eq.  (7),  define  M  =  [(?  0/[(l-  rR)nJ^].  For  a  wide  range  of  sky  and  sea  surface  conditions  and 
for  viewing  directions  relevant  to  remote  sensing,  JV/lies  in  the  range  of  0.53  to  0.55  [Mobley,  1994; 
Hoge  and  Lyon,  1996;  Hoge  et  al.,  1999a,b;  Morel  and  Gentili,  1996].  Thus  Mem  be  approximated 
as  A/  ~  0.54,  with  an  error  of  less  than  two  per  cent.  Using  Eq.  (7)  and  partitioning  the  beam 
attenuation  coefficient  as  c(z)  =  a{z)  +  b^z)  +  b^iz),  Eq.  (2)  becomes 


R 


rs 


Eda 

_ M/,(e,(p,0)  6,(0)/[27i  p/0)] _ 

-^0,(p,O)cos0  +  6^0)  [1  -4(0,9, 0)]  +  «(0)  +  V) 


(8) 


Except  for  the  small  error  associated  with  the  assumed  value  for  M,  Eq.  (8)  remains  an  exact  RTE 
expression  for  the  in-air  remote  sensing  reflectance,  R^,  just  above  the  sea  surface.  The  remote¬ 
sensing  reflectance  is  the  quantity  used  as  the  basis  for  ocean  color  remote  sensing  by  the  SeaWiFS 
[O’Reilly  et  al,  1998]  and  PRILLS  [Davis  et  al,  2002]  systems.  R,^  can  be  obtained  fi-om  at-sensor 
radiances  after  atmospheric  correction;  for  our  piuposes  here  it  is  therefore  considered  known.  As 
noted  by  Zaneveld  [1995],  it  is  desirable  to  use  RSRg,  the  in-air  value  of  RSR,  rather  than  R^  because 
the  scalar  irradiance  is  less  sensitive  to  solar  zenith  angle  effects  than  is  the  plane  irradiance  E^. 
We  thus  use  RSR^  =  R,,  to  rewrite  Eq.  (8)  as 


RSR  = 


"^oda 


M  (0) 

271  p^O) 


-^(0,9,O)cos0  +  ^>y(0)[l  -4(0,9,O)]  +  a(0)  +  6^(0) 


(9) 


The  simplicity  of  Zaneveld's  [1995]  original  in- water  RST?  formulation  remains  in  this  equation  for 
RSR^,  except  for  M  and  the  Pda/Pd  ratio  for  the  downwelling  light  field.  As  a  practical  matter,  Eq. 

(8)  is  presently  more  easily  applied  to  field  data  because  the  in-air  downwelling  plane  irradiance  is 
more  generally  available,  but  there  are  no  instrumental  barriers  to  using  scalar  irradiance  as  in  Eq. 

(9) . 
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To  further  simplify  Eq.  (9)  for  later  use,  define  the  first  term  in  the  denominator  as 


£>^(0,9,0)  =  -^0,9,0)  cosG.  (10) 

We  call  Dl  the  radiance  derivative  term  because  it  is  a  measure  of  the  depth  rate  of  change  of  the 
upwelling  radiance,  as  seen  in  Eq.  (5).  Define  the  second  term  in  the  denominator  as 

5/0,9, 0)  ^  6/0)  [1 -//0,9,O)] .  (11) 


The  shape  factory^  varies  from  0.963  to  1.152  for  nadir  viewing  [Weidemann  et  al.,  1995;  Zaneveld 
1995];  thus  5^  ranges  from  0.0376^to  -0.1526^  which  is  a  small  fraction  of  the  forward  scattering 
coefficient  bf.  We  therefore  call  5y^the  fractional forward  scattering  coefficient.  Thus,  bfmdfr  are 
found  in  a  combination  in  which  one  (/i)  serves  to  reduce  the  size  of  the  other  {bfr  Finally,  define 
the  mean  cosine  ratio  as 


R  = 

"  MO) 


(12) 


Using  definitions  (10)-(12),  Eq.  (9)  becomes 


M 


555/0^,9)  = 


//6.<P»0) 

2jt 


^6/0) 


^oda 


d/0,9,0)  +  5/0,9,O)  +  a(0)  +  6/0) 


(13) 


Equations  (8),  (9),  and  (13)  are  each  called  the  shape-factor  form  of  the  RTE;  we  shall  work  with 
Eq.  (13).  Our  ultimate  goal  is  to  use  Eq.  (13)  to  relate  the  unknown  absorption  and  backscatter 
coefficients  just  beneath  the  sea  surface  to  the  known  remotely  sensed  reflectance  and  other  known 
quantities.  As  noted  above,  M  ~  0.54.  However,  the  four  quantities  f,  R^,  D^,  and  [or, 
equivalently,^,  5^,  k,  and  6/(1-^)  as  seen  in  Eq.  (9);  for  brevity  we  refer  to  all  of  these  quantities 
as  shape  factors]  are  unknown.  The  shape  factors  depend  in  complicated  ways  on  the  water-column 
lOPs,  environmental  conditions  (sky  radiance,  sea  state),  and  viewing  geometry  (sun  zenith  angle, 
viewing  direction).  In  Section  4  we  show  how  to  model  the  shape  factors,  so  that  they  too  can  be 
considered  known  in  Eq.  (13). 

Note  finally  that  if  and  5^  are  both  zero,  then  both  and  RSR^  are  proportional  to  b^/ia  +  6b), 
which  is  a  well  known  approximate  functional  form  for  the  dependence  of  these  reflectances  on  the 
absorption  and  backscatter  coefficients. 
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3.  Linear  Form  of  the  Radiative  Transfer  Equation  and  Its  Inversion 


The  in-air  RSRj  of  Eq.  (13)  immediately  yields  the  fundamental  linear  form  of  the  RTE, 

a(0)  +  ^>*(0)  V  +  0,  (14) 


where 


F(0,9,O)  =  1  - 


2%  ^ 


(15) 


''oda 


V  is  called  the  backscattering  enhancement  factor  because,  in  general,  |  F|  >  10. 


We  next  partition  the  total  absorption  coefficient  into  contributions  by  pure  water, 
phytoplankton,  and  CDOM  plus  detritus.  Similarly,  the  backscatter  coefficient  is  written  as  the  sum 
of  contributions  by  pure  sea  water  and  by  particulate  matter.  It  is  easy  to  show  [Hoge  and  Lyon 
1996;  Hoge  et  al.,  1999a ;  Hoge  et  al.,  1999b]  that  the  equation  describing  the  desired  phytoplankton 
absorption  coefficient  CDOM+detritus  absorption  coefficient  aj,  and  total  constituent 
backscattering  coefficient  b^,  resulting  from  Eq.  (14)  is 

*  "A)  *  W  y  =  -oA)  -  r-D,-  s,.  (16) 


The  wavelength  dependency  of  the  lOP's  is  now  shown  explicitly,  while  the  depth  and  angular 
dependencies  have  been  suppressed  for  clarity.  Note  that  the  observed  water-leaving  radiances  Z,„„ 
occur  on  both  sides  of  the  equation  (within  V).  The  pure-water  absorption  a„  is  known  from  Pope 
and  Fry  [1997],  and  the  water  backscattering  coefficient  b^y/  is  given  by  Smith  and  Baker  [1981]. 
The  right  hand  side  of  Eq.  (16)  is  therefore  known,  given  the  shape  factors  and  a  measurement  of 
RSRJf^.  This  linear  form  of  the  RTE  is  still  exact  in  the  sense  that  no  approximations  have  been 
made  to  the  RTE,  but  clearly  the  lOP  retrieval  accuracy  will  be  determined  by  the  accuracy  of  the 
shape  factor  models. 

Given  the  water-leaving  radiance  at  three  wavelengths,  Eq.  (16)  still  cannot  be  solved  for  the  “big 
three”  lOPs,  apfkf  aff^  and  b^,(}^,  because  each  measurement  of  ^^^^(A,)  yields  an  equation  with 
three  unknown  lOPs.  However,  it  is  easy  to  show  that  a  consistent  solution  is  available  by 
introducing  spectral  models  for  apff),  afk)  and  b*,(^,)  [Hoge  and  Lyon  1996;  Hoge  et  al.,  1999a, 
Hoge  et  al.,  1999b].  Substitution  of  such  spectral  models  for  apfkf  aff)  and  b^^^  into  Eq.  (16) 


6 


yields 


ciph{\)^^V 


i\-Kf 


+  a^X^)exp[-5(X.-A,^)]  +  VV 


(  .  \ 


yh, 


V  = 


(17) 


-"A)  -  w  y-D,-  B,. 


Equation  (17)  now  has  only  three  unknowns,  ap^(Xg),  aj(Xj)  and  b^,(K)’  so  that  the  system  is  solvable 
given  measurements  of  RSRX^^i)  at  three  wavelengths.  This  linear  form  of  the  radiative  transfer 
equation  remains  exact  and  precise  but,  again,  the  uncertainty  in  the  retrieved  lOPs  is  now 
additionally  influenced  by  the  uncertainty  in  the  lOP  models  [Hoge  et  al.,  1999a]. 

At  their  respective  reference  wavelengths  Xg,  Xj ,  Aj,  the  lOPs  aJ7^)  and  linearly 

related  to  the  column  matrix,  or  vector,  containing  the  hydrospheric  constants  (sea  water  absorption 
and  backscattering),  radiances,  and  the  shape  factors.  [It  is  easy  to  see  from  Eq.  (16)  that  it  is  in 
principle  possible  to  concurrently  solve  for  the  radiance  derivative  term  Z)l(^)  and/or  for  the 
fractional  forward  scattering  coeflBcient,  in  addition  to  the  lOPs,  given  measurements  of  RSRJ}.^ 
at  additional  wavelengths.  However,  the  lOP  models  then  must  be  of  sufficient  accuracy  at  yet  a 
foiirth  and/or  fifth  wavelength,  and  the  required  wavelength  dependency  of  these  models  was  not 
a  focus  of  this  study.  Therefore,  such  retrievals  are  beyond  the  scope  of  this  initial  RTE  inversion 
work.]  Equation  (17)  is  very  similar  to  the  one  used  to  analyze  the  effect  of  radiance  errors  and 
model  uncertainties  upon  lOPs  [Hoge  and  Lyon,  1996],  and  to  retrieve  lOPs  from  airborne  upwelling 
radiances  [Hoge  et  al.,  1999a,  b]  when  retrieved  by  semianalytic  radiance  model  inversion.  Equation 
17  can  therefore  be  written  in  matrix  form  as 

Dp  =  h.  (18) 


Here  the  hydrospheric  vector  h  is  given  by  the  right  hand  side  of  Eq.  (17).  The  lOP  state  vector  is 
p  =  [aph(k^,  acfh),  h,Q%)V^  where  T  denotes  the  transpose,  and  D  is  the  data-model  matrix  [Hoge 
and  Lyon  1996;  Hoge  et  al.,  1999a,  b],  which  also  contains  shape  factors.  The  lOPs  are  immediately 
determined  from  p  =  D'  h.  (A  least  squares  solution,/?  =  [D^  /)]''  h,  can  be  used  if  the  number 
of  sensor  bands  exceeds  the  number  of  imknown  lOPs.) 

The  uncertainties  in  the  lOP  state  vector  p  can  be  analyzed  in  a  manner  similar  to  other  linear 
inversions  [Hoge  and  Lyon,  1996].  Since p  =  h,  both  D  and  h  determine  p  and  the  errors  that 
propagate  into  p.  Because  the  backscattering  shape  factor is  always  found  within  D,  influences 
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the  propagation  of  errors  into  the  lOPs  more  so  than  do  the  remaining  factors.  The  discussion  of  the 
uncertainties  in  the  lOP  state  vector  p  caused  by  possible  singularity  of  2)"'  and  by  perturbations  in 
D  somewhat  parallels  a  similar  previous  discussion  [Hoge  and  Lyon  1996]  and  is  therefore  provided 
in  Appendix  A. 


4.  Models  for  Shape  Factors  and  Related  Quantities 

Inversion  of  Eq.  (17)  using  remotely  sensed  ocean  color  data  requires  knowing  V,  D^,  and  Bj. 
These  quantities  in  turn  depend  on  the  shape  factors,^  and^^j  the  diffuse  attenuation  coefficient  k, 
and  the  mean  cosine  ratio  as  defined  above.  For  ease  of  comparison  with  previous  work  on  shape 
factors  [Weidemann  et  al.,  1995]  and  to  make  the  underlying  physics  as  transparent  as  possible,  we 
explicitly  modeiyj,^,  k,  and  rather  than  and  Bf.  For  notational  convenience,  we  let  X„  i  = 
1,2,3,  and  4,  denote^,^,  k,  and  respectively. 

Because  shape  factors,  diffuse  attenuation  functions,  and  mean  cosines  all  depend  on  the  ambient 
radiance,  they  depend  implicitly  on  the  solar  zenith  angle  and  viewing  direction,  as  well  as  on  the 
lOPs.  The  solar  angle  and  viewing  direction  are  known  in  any  particular  remote  sensing  situation; 
these  geometric  quantities  are  thus  available  for  modeling  the  Xj  in  terms  of  known  quantities. 
However,  the  lOPs  are  unknown.  If  we  restrict  ourselves  to  an  explicit  inversion  of  the  RTE  to 
obtain  the  lOPs,  we  must  exclude  the  lOPs  from  the  models  for  the  X^.  However,  if  we  perform  an 
implicit,  or  iterative,  inversion  of  the  RTE,  we  can  include  the  retrieved  lOPs  in  the  Xj  models,  for 
the  following  reason.  In  an  iterative  inversion,  we  start  with  an  initial  guess  for  the  Xj,  derived  either 
fi-om  models  that  do  not  include  lOPs  or  from  physical  intuition.  (For  example,  a  reasonable  initial 
guess  for  fi,  would  be  one,  the  value  corresponding  to  a  constant  phase  function.  Similarly,  =  1, 
k=0,  and  =  1  would  be  acceptable  initial  guesses.)  Using  the  initial  guesses  for  the  Xj,  the  RTE 
is  inverted  to  obtain  initial  values  for  ap,,(Xg),  a^h),  and  from  which  the  lOPs  a  =  a,  +  a„  = 

(fph  +  «</  +  and  bf,  =  bj,,  +  bj,^  can  be  obtained  at  all  wavelengths  via  the  lOP  models  seen  in  Eq. 
(17).  The  Xj  models  developed  below  are  based  on  an  assumed  phase  function  for  particle  scattering. 
Taking  the  particle  phase  function  as  known,  the  total  constituent  (particle)  backscatter  fraction  B, 
=  bjb,  is  also  known.  Thus  the  total  constituent  forward  scatter  coefficient  b^  can  be  obtained  fi-om 
the  recovered  6*,  and  B,:  b^  =  b,„{\  -  B,).  The  total  beam  attenuation  coefficient  is  then  known  from 
c  =  a  +  bjj  +  bi„  +  b„.  We  therefore  allow  ourselves  to  develop  models  for  the  Xj  that  depend  both 
on  the  known  geometrical  (viewing  direction,  solar  direction)  and  physical  (wind  speed,  wavelength) 
parameters,  as  well  as  on  certain  lOPs  (namely  a,  c,  and  6*). 


8 


4.1.  The  Database 


To  begin  our  analysis,  120  Hydrolight  [Mobley  and  Sundman,  2001a,  b]  runs  were  made  using 
its  lOP  model  for  Case  1  waters  and  the  Petzold  “average  particle”  phase  function  [Mobley  et  al., 
1993]  for  scattering  by  the  particles.  This  Case  1  lOP  model  is  a  two-component  model:  pure  water 
plus  “everything  else.”  The  non-water  absorption  and  scattering  coefficients  are  parameterized  in 
terms  of  the  chlorophyll  concentration  according  to  commonly  used  models  seen  in  Mobley  [1994, 
Eqs.  (3.27)  and  (3.40).]  The  input  for  these  runs  covered  a  wide  range  of  chlorophyll  concentrations, 
solar  zenith  angles,  cloud  covers,  and  wind  speeds.  Each  Hydrolight  run  generated  ouQiut  at  various 
wavelengths,  depths,  and  viewing  directions.  The  resulting  database  potentially  contains  millions 
of  records,  where  one  record  corresponds  to  a  particular  set  of  input  values,  output  values  for  a 
particular  viewing  geometry,  wavelength,  depth,  etc.,  and  the  values  of  the  four  X^.  Some  of  these 
records  are  not  of  great  interest,  e.g.  records  whose  azimuthal  viewing  directions  (p^  differ  by  only 
15  degrees  (the  resolution  of  in  the  standard  version  of  Hydrolight).  Therefore,  selected  records 
were  used  to  generate  a  database  of  more  manageable  size,  but  one  that  still  covers  the  range  of 
parameter  values  relevant  to  most  remote  sensing.  Table  1  shows  the  input  and  output  values  in  this 
database,  which  was  used  in  the  initial  investigation  of  the  functional  forms  of  the  Xj.  Each  of  the 
four  was  computed  for  each  parameter  combination  represented  in  Table  1. 


parameter 

values  used  in  Hydrolight  runs 

chlorophyll  concentration,  Chi 

solar  zenith  angle,  Oj 

0.0, 10, 20, 30, 40, 60  degrees 

wind  speed 

0, 10  m  s'' 

cloud  cover 

0, 100%,  i.e.,  clear  sky  and  solid  overcast 

wavelength,  X 

412, 426, 440, 465, 490, 522.5,  555,  612, 670, 
685  nm 

polar  viewing  angle  0„  (in  water,  relative  to 
the  zenith) 

0.0, 10, 20, 30, 40, 50, 60  degrees 

azimuthal  viewing  angle  cp^  (relative  to  sun) 

0, 90, 180  degrees 

depth,  z 

0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5  m 

Table  1.  Parameters  and  their  values  used  to  generate  the  original  database  of  184,800  records.  The 
values  shown  in  bold  correspond  to  the  1,500  records  in  the  remote-sensing  database. 
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We  first  examined  the  sensitivity  of  the  Xj  to  the  various  parameters  (wind  speed,  viewing 
direction,  lOPs,  etc)  available  for  constmction  of  models  for  the  X^.  We  found  that  the  surface  wind 
speed  has  a  negligible  effect  on  each  of  the  Xj  (less  than  one  percent  difference  in  Xf  for  the  0  and 
10  m  s  '  wind  speeds,  with  all  else  being  equal).  This  means  that  we  could  eliminate  the  wind  speed 
from  consideration  in  the  subsequent  modeling.  Likewise,  we  considered  only  the  clear-sky  data, 
which  are  of  greatest  interest  for  remote-sensing  applications.  As  noted  above,  for  remote  sensing 
applications  we  need  to  evaluate  the  shape  factors  Xj  only  at  depth  z  =  0.  (Note,  however,  that  the 
values  of  the  A!;  at  z  =  0  incorporate  the  effects  of  all  the  absorption  and  multiple  scattering  occurring 
throughout  the  entire  water  column.)  Thus  we  need  retain  only  the  output  from  the  Hydrolight  runs 
at  z  =  0.  The  original  database  included  records  generated  for  azimuthal  viewing  angles  of  (p^  =  0 
(looking  toward  the  sun)  and  (p^  =  180  degrees  (looking  away  from  the  sim).  Most  remote  sensing 
is  done  at  azimuthal  angles  of  (p^  «  90  degrees,  which  minimizes  sun  glint  and  instrument  self¬ 
shading.  Thus  we  retained  only  the  records  corresponding  to  (p^  =  90  degrees.  Likewise,  remote 
sensing  generally  uses  in-air  nadir  viewing  directions  Oyj  of  less  than  60  degrees,  which  corresponds 
to  in-water  angles  0^  <  40  degrees.  Eliminating  the  larger  in-water  nadir  viewing  angles  (0^  =  50  and 
60  degrees  in  Table  1)  gives  a  final  “remote  sensing”  data  set  of  1500  records,  which  was  used  to 
determine  models  for  the  Xj.  The  parameter  values  corresponding  to  this  remote-sensing  database 
are  shown  in  bold  in  Table  1. 


4.2.  Determining  Functional  Forms 

LetP^,  denote  the  parameters  (wavelength,  viewing  direction,  lOPs,  etc)  to  be  used 

in  modeling  the  These  parameters  include  those  seen  in  Table  1,  as  well  as  the  absorption, 
scattering,  and  backscattering  coefficients  (which  are  functions  of  the  chlorophyll  concentration  if 
Case  1  water  is  assumed). 


The  simplest  possible  model  for  the  Xj  is  a  linear  function  of  the  P*: 


N, 


X, 


-T, 

k=\ 


i  =  1,...,4. 


(19) 


The  are  fitting  coefficients  whose  values  are  to  be  determined;  a  different  set  of  coefficients  is 
needed  for  each  factor  A).  We  first  set  up  a  large  linear  least  squares  problem  to  see  if  this  model  was 
adequate  to  fit  the  various  factors.  Not  surprisingly,  the  fits  were  unsatisfactory.  In  other  words,  the 
world  is  more  complicated  than  Eq.  (19). 
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The  goal  then  became  to  develop  models  which  could  be  nonlinear  in  the  parameters,  but  which 
still  reflect  the  underlying  physics  of  the  Xj.  Thus  we  replace  Eq.  (19)  by 

m) 

=  E^,»(TO.K»-  (20) 

*=l 

This  notation  means  that  each  function  will  contain  some  subset  of  the  parameters  and  will 
have  its  own  set  of  fitting  coefficients.  For  example,  as  we  shall  see  in  the  next  section,  the  model 
for  fi,  will  contain  terms  involving  the  scattering-to-backscattering  ratio  and  a  nonlinear  function  of 
the  three  geometric  parameters  (0s,e^,(pj,  with  four  fitting  coefficients  in  all. 

4.3.  The  Model  for 

We  now  outline  the  development  of  the  model  for  To  obtain  initial  guidance  about  the 
possible  functional  form  for  an^^  model,  we  first  used  single-scattering  theory  to  evaluate  the 
definition  of/*  seen  in  Eq.  (3).  According  to  the  single  scattering  approximation  (SSA),  the 
downwelling  radiance  is  [Gordon,  1994,  Eq.  (1.30)] 

8(ip'  -  (p,)  e  ,  (21) 


where  p  =  cosO  and  (Ps>  9s)  denotes  the  in-water  direction  of  the  sun’s  direct  beam  as  refracted 
through  a  level  sea  surface;  6  is  the  Dirac  delta  function.  Inserting  (21)  into  (3),  integrating,  and 
evaluating  the  result  at  z  =  0,  gives 

4  27t  A  p(p^,(p^,p,(p) .  (22) 

"b 

Note  that  for  isotropic  scattering,  P  =  l/47t,  b  =  26*,  and y*  =  1,  as  expected.  As  the  phase  function 
becomes  more  anisotropic,  6/6*  increases,  and ^  increases.  In  Eq.  (22),  the  phase  function  tells  us 
that,  to  first  order,  y*  involves  downwelling  radiance  that  is  singly  scattered  from  the  sim’s  direct 
beam  into  the  viewing  direction.  In  remote  sensing  applications,  the  total  phase  function  (water  plus 
particles)  will  be  unknown,  but  for  any  given  phase  function  the  contribution  tof^  will  depend  on  the 
scattering  angle  9  corresponding  to  the  sun’s  downward  beam  being  scattered  into  the  upward 
viewing  direction.  This  scattering  angle  is  equivalent  to  the  easily  computed  sun-sensor  included 
angle  as  shown  in  Fig.  1.  Given  the  sun's  location  (0^,  (p^  =  0)  and  the  viewing  direction  (0^,  (pj, 
^  is  given  by 
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cos^  =  cos0^cos0y  +  sin05  sin0y  cos((Py) . 


Thus,  as  a  preliminary  functional  form  for  the ^  model,  we  consider 


where  S(^)  is  a  function  of  angle  ^  whose  form  is  to  be  determined. 

Figure  2  shows  the  1,500  values  of /*  in  the  remote  sensing  database  plotted  as  a  function  of  ^ 
and  blb^.  This  figure  suggests  that  a  cosine  function  may  capture  the  §  dependence.  Thus  we  try  a 
model  of  the  form 


/*  =  “i  a2  — cos(a3^). 


where  a„  Oj,  and  Uj  are  fitting  coefficients  (the  %  of  Eq.  (20)  for  model  /  =1)  whose  values  are  to 
be  determined  by  minimizing  the  squared  difference  between  and^*  for  the  1,500  values  in  the 
remote-sensing  database.  We  do  note  fi-om  the  color-coded  points  in  Fig.  2  that^  is  largest  for  small 
bib,,,  and  vice  versa,  which  is  contrary  to  the  behavior  predicted  by  Eq.  (22).  This  reversal  may  be 
due  to  the  dominance  of  multiple  scattering  in  ocean  waters,  but  further  investigations  would  be 
necessary  to  understand  this  discrepancy  between  the  SSA  predictions  and  the  Hydrolight 
predictions,  which  include  all  orders  of  multiple  scattering  and  other  effects  not  included  in  the  SSA. 
In  any  case,  there  is  a  clear  dependence  on  bib,,,  which  can  be  modeled. 
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The  best-fit  coefficients  aj  in  Eq.  (23)  were  determined  by  least-squares  minimization  using  a 
variety  of  numerical  techniques  appropriate  for  non-linear  functions.  After  comparing  the  model 
predictions  with  the  actuaiyj  values,  is  was  seen  that  not  all  of  the  b/b^  dependence  was  captured 
by  the  model  of  Eq.  (23).  Some  experimentation  showed  that  the  remaining  b/b^  dependence  could 
be  accounted  for  by  adding  another  term  proportional  to  b/b,,.  Thus  the  final^  took  the  form 


A  =  «i  a2-^cos(a3^)  +04-^ 


=  a,  +04 


a~ 

1  +  — cos(a3  4) 


(24) 


The  second  form  of  Eq.  (24)  shows  that  the  additive  term  is  equivalent  to  keeping  the  general  form 
of  the  model  suggested  by  the  SSA,  but  picking  a  different  angular  fimction  S(4).  The  complicated 
dependence  off^  on  b/b^  and  ^  is  not  surprising  if  one  remembers  that  the  S(^)  function  is 
fundamentally  an  attempt  to  parameterize  the  unknown  phase  fimction  effects  for  a  given  srm  and 
viewing  geometry;  thus  the  scattering  and  geometric  effects  are  not  independent.  The  final  set  of 
fitting  coefficients  for  model  (24)  is  shown  in  Table  2. 


coefficient 

value 

tti 

1.2077 

02 

0.001977 

03 

3.3790 

O4 

-0.004863 

Table  2.  Best-fit  coefficients  for  the^  model  of  Eq.  (24). 

Figure  3  shows  the  model  and  actual^  values.  The  dashed  lines  are  the  5%  error  bovmds.  Using 
the  model  of  Eq.  (24),  96.3%  of  the  predicted  values  are  within  5%  of  the  correct  value;  the  linear 
correlation  coefficient  between  the  model  and  actual  points  is  r  =  0.955.  There  is  no  systematic 
dependence  on  b/b^  or  ^  of  the  model  discrepancies  seen  in  the  individual  points  of  Fig.  3.  It  would 
be  possible  to  continue  adding  ad  hoc  terms  in  other  variables  to  Eq.  (24),  and  perhaps  reduce  the 
model-data  discrepancy  even  more.  However,  such  a  process  is  likely  to  deviate  fi"om  physical 
foundations,  with  the  end  result  that  the  final  model  would  not  be  applicable  beyond  the  exact 
conditions  used  to  generate  the  present  remote-sensing  data  set.  For  this  initial  study,  it  is  best  to 
be  content  with  the  simple  model  of  Eq.  (24). 
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Fig.  3.  Comparison  of  modeled  and  actual  ^  values,  using  the  model  of  Eq.  (24).  96.3%  of  the 
modeled  values  lie  within  the  dashed  lines,  which  represent  values  with  5%  of  the  correct  value;  the 
model-actual  correlation  coefficient  is  r  =  0.955. 


4.4.  The  Model  for^ 

Just  as  with^,  we  first  tried  to  use  the  SSA  for  guidance  as  to  the  general  form  of  the  model  for 
fi.  In  the  SSA,  the  upwelling  radiance  is  [Gordon,  1994,  Eq.  (1.32)] 

L  inWyZ)  =  -^/0)p(p^,(p^,p',(p') — 1— e  (25) 

c  p,  -  p' 


Note  that  p^  >  0  and  p'  <  0.  Inserting  this  SSA  radiance  into  Eq.  (4),  the  definition  off^,  integrating, 
and  setting  z  =  0  gives 


^  f  f  P(p',(p',p,(p)  P(p,,(p^,p',(p') — — 

bfcV:°  {  {  M.-li 
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In  most  ocean  waters,  b  ~  bf.  However,  further  simplification  is  difficult.  The  remaining  integrals 
describe  how  the  sun’s  downwelling  direct  beam  is  first  scattered  upward,  and  then  scattered  again 
into  the  viewing  direction.  The  most  we  can  say  is  that  this  is  some  function  of  the  scattering  phase 
function  and  the  viewing  geometry.  As  with^,  we  tried  to  parameterize  this  function  in  terms  of  the 
sun-sensor  included  angle  ^  via  a  model  of  the  form 


Unfortunately,  plots  of^  similar  to  Fig.  2  did  not  suggest  a  clear  functional  form  for  H(4)  or  show 
any  significant  dependence  on  b/c.  This  failure  of  the  SSA  to  provide  a  functional  form  for^  is  not 
surprising  because  inherently  involves  at  least  two  scatterings,  and  multiple  scattering  can  be 
expected  to  make  an  important  contribution  to  the  upwelling  radiance. 

We  next  examined  the  linear  correlations  between^^  and  the  various  available  fitting  parameters. 
The  results  are  seen  in  Table  3. 


parameter 

correlation  coef.  r 

with/t 

absorption  coefficient  a 

0.291 

scattering  coefficient  b 

0.205 

albedo  of  single  scattering  b/c 

-0.076 

backscatter  coef  bj, 

0.180 

forward  scatter  coef  bf 

0.206 

backscatter  fraction  bf/b 

-0.343 

forward  scattering  fraction  b/b 

0.343 

backscatter  to  absorption  b^/a 

-0.231 

wavelength  X 

0.210 

solar  zenith  angle  Gj 

0.778 

polar  viewing  angle  0^ 

-0.006 

sun-sensor  included  angle  ^ 

-0.562 

Table  3.  Correlation  between  the  radiance  shape  factor^  and  various  parameters. 
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The  only  potential  model  parameter  that  correlates  even  moderately  strongly  with^i  is  the  solar 
zenith  angle  G^.  A  plot  of/^  vs.  0^  suggested  a  sine  function  for  0,  (although  a  simple  linear  function 
is  almost  as  good).  Thus  we  first  tried  a  model  of  the  form  =  a,  +  Oj  sin(a30^) .  The  residuals 
of  this  model  showed  a  weak  wavelength  dependence.  After  considerable  experimentation,  we 
settled  on  the  model 


(26) 


The  best-fit  values  of  the  coefficients  are  a,  =  1.0247,  Oj  =  0.4584  (for  X,  in  nanometers),  and  03= 
0. 1221  (for  0s  measured  in  radians).  Figure  4  shows  the  scatter  plot  for  this  model.  93.0%  of  the 
model  predictions  are  within  2%  of  correct  (points  lying  between  the  dashed  lines);  the  correlation 
coefficient  is  r  =  0.829.  Although  it  is  possible  to  obtain  slightly  better  fits  by  including  lOPs  in  an 
ad  hoc  fashion,  we  choose  to  stop  with  the  model  of  Eq.  (26)  because  of  its  simplicity  and  because 
of  the  lack  of  physical  guidance  for  the  lOP  dependence. 


Fig.  4.  Comparison  of  modeled  and  actual  values,  using  the  model  of  Eq.  (24).  93.0%  of  the 
modeled  values  lie  within  the  dashed  lines,  which  represent  values  with  2%  of  the  correct  value;  the 
model-actual  correlation  coefficient  is  r  =  0.829. 
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Although  we  were  unable  to  model  the  remaining  variability  offi  in  terms  of  the  lOPs  or  other 
parameters,  this  may  be  of  little  importance  in  predicting/^  itself,  because/  is  always  near  one. 
However,  the  variability  in  /  determines  the  variability  in  the  fractional  forward  scattering 
coefficient  Bf=  Zy(l  -/).  Small  firactional  errors  in/  can  cause  large  fi^ctional  errors  in  Bp  Figure 
5  shows  the  resulting  scatter  plot  for  Bp  computed  using  the  exact  values  of  bj  as  found  in  the 
database.  Although  93.0%  of  the/  values  are  within  2%  of  their  correct  value,  only  5.4%  of  the  B^ 
values  are  within  2%  of  the  correct  value;  58. 1%  of  the  Bfdxt  within  20%  of  correct.  However,  a 
percentage  error  criterion  may  be  misleading  for  because  of  the  cluster  of  points  near  zero,  where 
small  absolute  errors  can  be  large  firactional  errors.  The  dotted  lines  in  Fig.  5  thus  show  absolute 
errors  of  ±0.02  m'*;  94.5%  of  the  Jyhave  errors  smaller  than  this.  The  model-actual  correlation 
coefficient  is  r  =  0.966. 


Fig.  5.  Comparison  of  modeled  and  actual  values,  computed  using  the/  model  of  Eq.  (26).  The 
dashed  lines  are  20%  error  bounds;  the  dotted  lines  are  ±0.02  m  '  error  bounds. 
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4.5.  The  Model  for  ^ 


The  SSA  again  suggests  a  functional  form  for  the  ^0,(p,z=O)  model.  Differentiating  the  SSA 
upwelling  radiance  of  Eq.  (25)  gives 


k(Q,(p,0) 


1 

L„(z)  dz 


c 

cos(0^) 


Although  k  is  strongly  correlated  with  a{r  =  0.99),  it  is  not  strongly  correlated  with  c(r  =  0.56).  At 
the  level  of  the  quasi-single-scattering  approximation  (QSS  A),  which  often  works  well  for  upwelling 
radiances,  c  «  a  +  h*,  in  which  case  A: »  (a  +  hj)/cos(0s).  This  suggested  that  we  first  try  to  model  k 
with  the  functional  form 


^  ^  «i  ^  +  a2 
cos(0^) 


(27) 


When  Eq.  (27)  was  used  to  fit  the  points  in  the  remote  sensing  data  set,  the  points  separated  into 
distinct  groups  for  0^  ^  60°  and  for  0^  <  60°.  We  then  plotted  k  and  the  residual  error  in  k  as 
functions  of  0^,  which  suggested  that  an  additive  sin(0s)  term  would  represent  the  0^  dependence 
better  than  the  l/cos(0J  factor  in  Eq.  (27).  This  then  gave  the  final  k  model: 

^  =  a,  a  +  ttj  +  ttj  sin(0^) .  (28) 

The  best-fit  parameters  are  Uj  =  1.0896,  -0.5931,  and  Uj  =  0.0492  (for  a,  b/,,  and  k  in  inverse 

meters).  [Note  that  although  0,  in  Eq.  (28)  refers  to  the  polar  angle  of  the  solar  beam  in  water,  we 
can  use  the  in-air  solar  zenith  angle  for  0^,  because  the  index-of-refiraction  factor  that  converts  sin(0s 
in  air)  to  sin(0s  in  water)  is  incorporated  into  a^]  The  points  still  separate  somewhat  by  0^  ^  60°  and 
for  0s  <  60°,  but  not  as  much  as  for  Eq.  (27).  Although  Eq.  (28)  has  lost  some  of  its  intuitive,  first- 
order  physics,  namely  the  l/cos(0s)  factor  in  Eq.  (27),  the  final  model  does  a  better  job  of  predicting 
k,  which  of  course  is  influenced  by  multiple  scattering  and  other  effects  not  included  in  the  SSA. 
Figure  6  shows  the  scatter  plot  for  the  k  model.  Because  of  small  differences  near  A:  =  0,  only  76.8% 
of  the  points  are  within  20%  of  the  correct  value.  However,  93.6%  of  the  points  lie  within  an 
absolute  error  of  ±0.05  m"'.  The  model-actual  correlation  coefficient  is  r  =  0.993. 
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4.6.  The  Model  for 


A  first-order  model  for  =  |t^in  air)  /  water  at  2=0)  can  be  obtained  simply  by  using 
Snell’s  law  to  refract  the  direct  solar  beam  through  a  level  water  surface.  The  result  is 


=  «i 


COS0„ 


cos 


sm 


'  sinO^'" 


/J 


=  a,  0(0^), 


(29) 


where  n=  1.34  is  the  index  of  refraction  of  water.  The  ratio  R^^  is  independent  of  the  viewing 
direction  because  the  mean  cosines  are  computed  from  integrals  over  direction  of  the  radiance 
distribution;  there  are  consequently  many  fewer  distinct  points  in  the  data  set.  When  Eq.  (29)  is 
applied  to  the  RS  data  set,  the  best-fit  value  of  a,  is  0.869.  Figure  7  shows  the  results  of  this  model 
applied  to  all  points  in  the  RS  data  set.  The  six  groups  of  points  correspond  to  the  six  solar  angles 
in  the  data  set:  O^  =  0,  10,  20,  30,  40,  and  60  degrees.  85.7%  of  the  points  lie  within  5%  of  the 
correct  value,  as  shown  by  the  dashed  lines.  The  correlation  coefficient  is  r  =  0.964. 


)U  . . I . . 

0.50  0.60  0.70  0.80  0.90  1.00 


actual 

Fig.  7.  The  model  of  Eq.  (29)  applied  to  the  RS  data  set.  The  dashed  lines  are  5%  error  boimds. 
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5.  Discussion 


We  now  have  developed  models  that  predict  the  shape  factors  and  related  quantities  with  vaiying 
degrees  of  certainty.  The  question  next  arises  as  to  how  well  we  can  predict  the  remotely  sensed 
reflectance  RSR^  using  these  models  in  Eq.  (9).  Figure  8  gives  the  answer:  68.1%  of  the  RSR^ 
predictions  fall  within  10%  of  the  correct  values,  and  95.7%  fall  within  20%  of  correct  (shown  by 
the  dashed  lines  in  Fig.  8);  85.6%  of  the  predictions  are  within  ±0.0005  sr"'  of  the  correct  value.  The 
model-actual  correlation  coefficient  is  r  =  0.983. 

We  have  already  noted  that  inversion  of  the  shape  factor  form  of  the  RTE  is  exact  from  the 
standpoint  of  radiative  transfer  theory  and  that  the  uncertainties  in  the  retrieved  lOPs  are  due  only 
to  the  accuracy  of  the  lOP  models  [Hoge  and  Lyon  1996;  Hoge  et  al.  1999a;  Hoge  et  al.,  1999b]  and 
to  the  models  for  the  shape  factors  and  related  quantities.  The  lOP  models  are  fairly  mature 
compared  to  the  preliminary  shape  factor  models  presented  here. 

To  understand  how  the  remaining  errors  in  the  shape  factor  models,  and  in  our  ability  to  predict 
RSR^,  affects  the  accuracy  of  the  retrieved  lOPs,  it  is  necessary  to  study  performance  of  the  iterative 
inversion  algorithm  outlined  above  via  its  application  to  both  synthetic  and  actual  data  sets.  The 
overall  inversion  is  investigated  first  in  the  controlled  environment  of  Hydrolight-generated  synthetic 
RSR^  data,  for  which  the  correct  lOP  values  are  known  from  the  input  to  Hydrolight.  The  inversion 
is  then  applied  to  actual  RSR„  field  data,  for  which  the  field  measurements  of  the  retrieved  lOPs  are 
themselves  subject  to  experimental  errors.  Those  evaluations  of  the  overall  inversion  algorithm  are 
the  subject  of  the  companion  paper  [Hoge  et  al.,  2002]. 
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Appendix  A.  Uncertainties  in  the  lOP  State  Vector  p 


Sensitivity  of  p  to  Perturbations  in  D 

Perturbations  within  D  arise,  for  example,  from  the  water-leaving  radiances,  scalar  irradiances, 
backscattering  shape  factor,  and  from  the  lOP  models.  Similarly,  imcertainties  in  h  arise  from  the 
radiances,  irradiances,  hydrospheric  constants  (or  lOP  constants  and  for  sea  water,  as  well 
as^,  dLJk)/(k,  and  cos  d.  Relative  to  h,  the  data-model  matrix,  D,  plays  the  major  role 

in  the  propagation  of  errors  into  p  since  ||  p  -  p'll  /||p|l  ^  k(D)  (||  A||  /1|D||  +  ||  8||  /||h||)  where  k(D) 
=llDi|  IID'^II  [Ortega  1990;  Hoge  and  Lyon  1996].  The  latter  expression  is  the  condition  number  of 
D,  and  A  and  8  represent  uncertainty  or  perturbation  of  D  and  h  respectively,  p'  is  the  perturbed 
solution  of  p.  The  first  expression  simply  states  that  to  first  order  the  relative  error  in  p  can  be  k(D) 
times  the  relative  error  in  D  and  h.  Thus,  the  propagation  into  p  of  the  relative  errors  of  both  D  and 
h  is  governed  by  the  condition  number  of  D.  For  any  norm,  1  ^  k(D)  <  <».  For  the  limiting  cases: 
k(D)  =  1,  D  is  said  to  be  perfectly  conditioned  while  for  k(D)  =  oo,  D  is  singular.  For  intermediate 
values  of  k(D)  the  interpretation  of  the  condition  number  is  very  subjective  and  must  be  evaluated 
separately.  For  "large"  k(D)  the  D  matrix  is  said  to  be  ill-conditioned  and  large  errors  may  be  found 
in  p.  For  "small"  k(D)  the  D  matrix  is  said  to  be  well-conditioned  and  smaller  errors  may  be  foimd 
in  p.  Of  the  five  components,  only  occurs  in  D  (via  V)  and  therefore  provides  the  strongest 
influence  on  the  lOP  retrieval  errors.  This  is  in  agreement  with  Zaneveld  [1995]  who  concluded  that 
yj  is  most  critical  since  the  in-water  remotely  sensed  reflectance  (see  equation  1)  is  directly 
proportional  to  it.  It  is  for  this  reason  that  RTE  component  model  developments  should  probably 
focus  onf^.  If  other  RTE  components  such  as  (or  dLJdz)  and/or  Bf  are  solved-for ,  then  they  too 
will  appear  in  the  D  matrix  and  thereby  further  increase  the  errors  in  p.  In  part,  the  additional 
uncertainty  in  p  will  be  due  to  and/or  19^  model  errors.  Also,  in  general  the  condition  number 
increases  as  the  number  of  unknowns  increases  [McCormick  1992]  contributing  to  additional 
uncertainty  in  p. 

Sensitivity  of  p  to  Perturbations  in  h 

While  the  condition  of  D  is  most  important  in  determining  the  errors  in  the  lOP  state  vector  p, 
the  error  propagation  equation,  I  p  -  p'||  /||p||  ^  ||D||  |D''||  (||  A||  /||D||  +  ||  81  /||h||)  shows  that 
uncertainties  in  h  also  propagate  into  p. 
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Notation 


i 


^ph 

a, 

Cly^f 

b 

h 

bf 

h. 

^bw 

Bf 

c 

CDOM 

D 

detD 

D, 

Eji^) 

Eoda 

E/z) 

Ed. 

E. .(z) 

Fik 

fb 

h 

A 

A 

fi 

fi 

g 

h 

lOP 


total  absorption  coefficient,  a,  +a„,  (m'');  denotes  “in  air”  when  used  as  a  subscript 

absorption  coefficient  of  CDOM  and  detritus  (m'*) 

absorption  coefficient  of  phytoplankton  (m'') 

total  constituent  absorption  coefficient,  a,  =  (m"') 

absorption  coefficient  of  water  (m  ') 

total  scattering  coefficient  (m'') 

total  backscattering  coefficient,  b,,  =  b^„  +  b^,  (m'') 

total  forward  scattering  coefficient  (m  ') 

total  constituent  backscattering  (TCB)  coefficient  (m"') 

backscattering  coefficient  of  seawater  (m'') 

fractional  forward  scattering  coefficient,  defined  by  Eq.  (9)  (m'‘) 

beam  attenuation  coefficient,  c  =  a  +  b  (m'‘) 

chromophoric  dissolved  organic  matter 

data  and  model  matrix 

determinant  of  D 

radiance  derivative  term,  defined  by  Eq.  (8)  (m'') 

downwelling  scalar  irradiance,  in  water  (W  m‘^  nm‘*) 

downwelling  scalar  irradiance,  in  air,  just  above  sea  surface,  (W  m'^  nm"') 

downwelling  plane  irradiance,  in  water  (W  m'^  nm"') 

downwelling  plane  irradiance,  in  air  (W  m‘^  nm"') 

upwelling  plane  irradiance,  in  water  (W  m'^  nm‘') 

upwelling  plane  irradiance,  in  air  (W  m'^) 

modeling  function 

backscattering  shape  factor,  dimensionless 

estimated  backscattering  shape  factor,  dimensionless 

radiance  shape  factor,  dimensionless 

estimated  radiance  shape  factor,  dimensionless. 

average  of  f^  values  having  0^  in  remote  sensing  data  set 

phytoplankton  Gaussian  model  spectral  width  parameter  (ran) 

vector  of  hydrospheric  constants,  shape  factors,  radiance  attenuation  coef.,  (m'*) 

inherent  optical  property 

diffuse  attenuation  coefficient  for  upwelling  radiance,  =  -d[logL„(r,0,(p)]/dz 
upwelling  radiance,  below  sea  surface  (W  m'^  sr'  ran’) 
upwelling  radiance,  in  air,  just  above  sea  surface  (W  m'^  sr '  ran  ’) 
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M 
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n„ 

Pk 

P 

r 

R 

R. 

RSR 

RSR^ 

Ks 

RTE 

S 

t 

r 

TCA 

TCB 

V 

z 

a, 

P 

Yi 

A 

4 

4 

4 

/4 

^da 

(P 

d 

% 

®o 


[(t  T)/[(l-rR)nJ]  ~  0.529  for  nominal  sea  conditions,  dimensionless. 

total  constituent  backscattering  coefficient  spectral  model  exponent,  as  used  in  Eq. 

(14),  dimensionless. 

index  of  refi:action  of  sea  water,  dimensionless, 
modeling  parameter 

oceanic  state  vector  of  retrieved  lOPs,  m  ’ 

water-to-air  surface  reflectance  for  plane  irradiance,  dimensionless 
plane  irradiance  reflectance,  in  water,  E,/z)/E/z),  dimensionless 
ratio  of  air-to-water  mean  cosines  for  downwelling  irradiance 
remotely  sensed  reflectance,  in-water,  L,/9,f,z)/Eo/z),  (sr') 
remotely  sensed  reflectance,  in  air,  L„/6,^)/E„j^  =  /?rs  (sr'‘) 
remote  sensing  reflectance,  in  air,  L,Jd,(p)/Ea^,  (sr'‘) 

Radiative  Transfer  Equation 

spectral  slope  within  the  model  for  CDOM/detritus,  (nm"') 
water-to-air  radiance  transmittance,  dimensionless 
water-to-air  plane  irradiance  transmittance,  dimensionless 
total  constituent  absorption,  a, ,  (m"') 
total  constituent  backscattering,  bbt ,  (m'') 

backscattering  enhancement  factor,  defined  by  Eq.  (12),  dimensionless 
shape  factors  and  related  quantities;  defined  in  Table  1  for  i  =  1,2,. ..6. 
depth,  (m) 

model  fitting  parameters,  see  tables  for  numerical  values 
volume  scattering  function 

model  fitting  parameters,  see  tables  for  numerical  values 
wavelength,  (nm) 

reference  X  for  total  constituent  backscattering  (TCB)  coefficient  model,  (nm) 

reference  wavelength  for  CDOM/detritus  absorption  coefficient  model,  (nm) 

peak  wavelength  for  Gaussian  phytoplankton  absorption  coefficient  model,  (nm) 

peak  wavelength  for  Gaussian  model  for  depth  derivative  of  radiance,  (nm) 

wavelength  of  observational  bands,  i=  1,  2,  3, ...,  (nm) 

average  cosine  for  downwelling  irradiance,  in  water,  dimensionless. 

average  cosine  for  downwelling  irradiance,  in  air,  dimensionless. 

azimuth  angle,  radians;  subscripts  v,  s  denote  solar,  viewing 

polar  zenith  angle  in-water,  radians;  subscripts  v,  s,  a  denotes  solar,  viewing,  in-air 

included  angle,  solar-to-viewing  direction 

single-scattering  albedo 
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